% 定义参数
omega = 1;  % 角频率

% 定义时间跨度
tspan = [0 1];

% 定义初始条件的范围和数量
x0_range = linspace(-5, 5, 50);  % 位置初始条件
v0_range = linspace(-5, 5, 50);  % 速度初始条件

% 创建图形和坐标轴
figure;
ax = axes;
hold on;
grid on;
xlabel('Position (x)');
ylabel('Velocity (dx/dt)');
title('Phase Space Plot of Simple Harmonic Oscillator');

% 定义微分方程
dxdt = @(t, x) [x(2); omega^2 * x(1)];

% 循环不同的初始条件
for x0 = x0_range
    for v0 = v0_range
        % 定义初始条件
        initial_conditions = [x0; v0];
        
        % 使用ode45求解微分方程
        [t, x] = ode45(dxdt, tspan, initial_conditions);
        
        % 绘制相空间图
        quiver(x(1,1), x(1,2), (x(end,1)-x(1,1))/5, (x(end,2)-x(1,2))/5, 0.5,'Color', rand(1,3));
    end
end

% 调整视图
axis equal;
xlim([-2 2]);
ylim([-2 2]);
hold off;